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ABSTRACT 

We jointly constrain the luminosity function (LP) and black hole mass function (BHMF) of broad-line 
quasars with forward Bayesian modeling in the quasar mass-luminosity plane, based on a homogeneous sam- 
ple of ~ 58,000 SDSS DR7 quasars at z ^ 0.3-5. We take into account the selection effect of the sample 
flux limit; more importantly, we deal with the statistical scatter between true BH masses and FWHM-based 
single-epoch virial mass estimates, as well as potential luminosity-dependent biases of these mass estimates. 
The LF is tightly constrained in the regime sampled by SDSS, and makes reasonable predictions when extrap- 
olated to ^ 3 magnitudes fainter Downsizing is seen in the model LF. On the other hand, we find it difficult to 
constrain the BHMF to within a factor of a few at z > 0.7 (with Mgii and Civ-based virial BH masses). This 
is mainly driven by the unknown luminosity-dependent bias of these mass estimators and its degeneracy with 
other model parameters, and secondly driven by the fact that SDSS quasars only sample the tip of the active BH 
population at high redshift. Nevertheless, the most likely models favor a positive luminosity-dependent bias for 
Mgii and possibly for Civ, such that at fixed true BH mass, objects with higher-than-average luminosities have 
over-estimated FWHM-based virial masses. There is tentative evidence that downsizing also manifests itself in 
the active BHMF, and the BH mass density in broad-line quasars contributes an insignificant amount to the total 
BH mass density at all times. Within our model uncertainties, we do not find a strong BH mass dependence 
of the mean Eddington ratio; but there is evidence that the mean Eddington ratio (at fixed BH mass) increases 
with redshift. 
Subject headings: black hole physics — galaxies: active — quasars: general — surveys 



1. INTRODUCTION 

One major effort in modern galaxy formation studies is 
to understand the cosmic evolution of supermassive black 
holes (SMBHs), given their ubiquitous existence in almost 
every local bulge-dominant galaxy, and possible roles during 
their co-evolution with the host gal axy (e.g., Magorrian et aQ 
19981: iGebhardt et alJ [2000; .Ferrarese & Merritt ,200a 
Gultekin et all 120091; iHopkins et all 120081; ISomervil le et alj 
2008!). Over the past decade, the rapidly growing body 
of observational data and numerical simulations have led 
to a coherent picture of the cosmic formation and evolu- 
tion of SMBHs withi n the hierarchical ACDM paradigm 
(e.g., [H aiman & LoeJ [19981 iKauffmann & Haehnelt 2000; 



'2002'; 'Yu & Lu' '200?, "loos'; 
Marconi etal. 2004; Merloni 



'Shankar et al.' '2004', 



Wvithe&Loeb 200* 'Volonteri et al. 2003; Hopkins et aj] 
2006, 2008; Shankar et al. 2009, 2010; Shen 2009)^ Although 
many fundamental issues regarding SMBH growth still 
remain unclear (such as BH seeds, fueling and feedback 
mechanisms), these cosmological SMBH models are starting 
to reproduce a variety of observed SMBH statistics in an 
unprecedented manner. 

It is now widely appreciated that SMBHs grow by gas 
accretion in the past, during which they are witn essed as 
quasars and active galactic nuclei (AGNs) (e.g., ISalpeted 
fl964tlZerdovich & Novikovll964l:li^den-BeUl969h . In the 
local Universe, the mass function of dormant SMBHs is esti- 
mated by convolving the galaxy distribution functions with 
various scaling relations between galaxy properties and BH 
mass. This relic SMBH population has been used to con- 
strain the accretion history of their active counterparts, us- 
ing the Soltan argument and its extensions (e . g.,[Sp ltan 1982; 
ISmall & Blandfordl[T992l:rSalucci et al.l[T999l: [Yr& Tremaind 
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2004; Hopkins etal. 2001 
iMerloni & Heinz 2008). The agreement between the relic BH 
mass density and the accreted mass density provides com- 
pelling evidence that these two populations are ultimately 
connected. Therefore it is of imminent importance to quan- 
tify the abundance of active SMBHs as a function of redshift. 

The demographics of the active SMBH population has been 
the central topi c for quasar s t udies since the first discov- 
ery of quasars (ISchmidtlll963L Il968h . Traditionally this is 
done in terms of the luminosity function (LF), i.e., the abun- 
dance of objects at different luminosities. Measuring the LF 
and its evolution has be en the most important goal for mod- 
ern qu asar surveys (e.g..l Schmidt & Greenlll983l ; lGreen et alj 
119861) . In the last decade the LF has been measured for differ- 
ent populations of active SMBHs and in different bands (e.g.. 
Fan etal. 2001, 2004; Bovle et al. 2000; Wolf etal. 2003t_ 
Croometal. 2004, 2009; ' Hao et al. 2005; Richard s et al 
, 2005 2006; Jia ng et al.. .2006. .2008. .2009; .Fqn tanot et al 
2007t [Bongiorn o et alj|2607t IWillott et al.ll2010bl; i Ueda et aT 
■2003^ Hasinger 'erarri2005l: ISiiverman et al.1 I2005l i2008r 
Barger et al. 2005|), and it constitutes a crucial observational 
component in all cosmological SMBH models. 

A more important physical quantity of SMBHs is BH mass. 
BH mass is directly related to growth, and when the BH is 
active, it determines the accretion efficiency (M/Mbh) via 
the Eddington ratio and an assumed radiative efficiency (e.g., 
such as the average value constrained by the Soltan argu- 
ment). Thus knowing the mass function of SMBHs as a func- 
tion of redshift adds significantly to our understandings of 
their cosmic evolution. 

It remains challenging to directly measure the dormant 
BHMF at high redshift. This is not only because the 
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galaxy distribution functions are less well-constrained at 
high redshift, but also because the evolution of the scal- 
ing relations (both the mean relation and the scatter) be- 
tween galaxy properties and BH mass is poorly understood. 
On the other hand, it has become possible to measure the 
active BHMF of broad-line quasars^, using the so-called 
virial BH mass estimators based on their broad emission 
line and continuum propert ies m e asured from si ngle-epoch 
spectra (e.g., |W andel et all 119991; iMcLure & Du nlop 2004; 
IVestergaar d & PetersonI 12006). a technique rooted on rever- 
beration mapping (RM) studies of local broad-line AGNs 
(e.g. JBlandford & McKee 1982; Peterson 1993; Kas pi et all 
I2000t [Peterson et al..i2004; iBentz et al.M2006t i2009a)." These 
single-epoch virial BH mass estimators are calibrated em- 
pirically using the RM AGN sample to yield on average 
consistent BH mass estimates compared with RM masses, 
which are further tied to the BH masses predicted using 
the M bh ~ g relation (e.g.. iTremaine et al]|2002l: lOnken et alj 
l2004l) . The nominal scatter between these single-epoch virial 
estimates and the RM mass es i s on the order of ;^ 0.4 
dex (e.g., McLure & Jarvis 2001 iMcLure & Dunlop I l2004t 
IVestergaard & Petersonli2006t 

A couple of recent studies have applied this technique to 
measure the a ctive BHMF with st a tistical quasar and AGN 
samples (e.g. Greene & Ho '2007"; Vestergaard et al. [20081: 
IVestergaard '& Osmer 2009; Schulze & Wisotzki 201(|. A 
robust determination of the active BHMF constitutes an im- 
portant building block of cosmological SMBH models, in 
addition to the luminosity function. These virial mass es- 
timators also enable statistical studies on the Eddington ra- 
tios of broad-line quasars and AGNs (e.g., Vestergaard 200^ 
McLure & Dunlop 2004; Kollmeier et al. 2006; Sulentic et aD 
2p06; Babic et al.. .2007.: Jiang et al. 2007; K urk et al.. ^OOtT 
Netzeretal. 20 07t IShen e t al. 2008; Gavig naudet alll2008[ 
Labita 2009; Trump et al.i i2009L i2011i: iWillott et al.ll2010ai 
Trakhtenbrot et al. 2011), over a wide range of luminosities 
and redshifts, and therefore provide constraints on the accre- 
tion efficiency of these active SMBHs. 

With the development of these virial mass estimators, we 
now have both BH mass estimates and luminosities for broad- 
line quasar samples. Given the intimate relation between BH 
mass and luminosity, it is important and necessary to study 
their joint distribution and evolution in the mass-luminosit y 
plane (e.g.. lSteinhardt & ElvisEolOatlSteinhardt et al.ll201 Ih . 
This represents a significant step forward to study the demog- 
raphy of quasars than using LF alone, and offers new insights 
on the properties and evolution of the active SMBH popula- 
tion. 

However, the importance of distinguishing between virial 
mass estimates and true BH masses can hardly be over- 
stressed. While these virial estimators currently are the 
only practical way to estimate BH masses for large sam- 
ples of broad-line quasars and AGNs, the nontrivial uncer- 
tainty of these imperfect estimators has severe impact on 
the mass distribution under study. The difference between 
virial masses and true masses not only modifies the underly- 
ing true distrib ution, but also introduce s Malmquist-type bi- 
ases (e.g., Shenetal.1 120081: iKellv et al.i 2009; Shen & Kellvl 
I2OIOI) . These effects tend to dilute an y potential mass- 
dependent trends or correlations (e.g., iKellv & Bechtoldl 
|2007; She n et al]l2009l) . and may lead to unreliable conclu- 

■^ From now on, unless otherwise specified, we use the term "quasar" to 
refer to broad-line (type 1 ) quasars for simplicity. 



sions. Thus it is important to consider these effects when the 
st atistics is be coming good enough. 

iKellv et al] (|2()09) developed a Bayesian framework to es- 
timate the BHMF/LF for broad-line quasars, which accounts 
for the uncertainty in virial BH mass estimates, as well as 
the selection incompleteness in BH mass (since the sample 
is selected in luminosity). This method was subsequentl y 
applied to the SDSS DR3 quasar sample (iRellv et al."2010!), 
based on virial mass estimates from Vestergaard et al. (2008). 
This Bayesian framework is a more rigorous and quantita- 
ti ve treatmen t than the simple forward modeling performed 
in lShen et al.l (l2008l) . and allows a more reliable measurement 
of the true active BHMF and its uncertainty for quasars. 

Equipped with an improved version of this Bayesian frame- 
work, in this paper we measure the active BHMF and LF 
based on a homogeneous sample of ~ 58,000 quasars from 
SPSS DR7 with FWHM-based virial mass estimates from 
IShen et aLJ ( 12011^ . The much improved statistics now allows 
a detailed examination of the joint distribution in the mass- 
luminosity plane, and provides better constraints on BH ac- 
cretion properties. 

A key difference in our approach compared with most ear- 
lier work is the attempt to account for the uncertainty (error) 
in these virial mass estimates. We distinguish three types of 
errors in single-epoch virial BH mass estimates: 

• measurement error, which is propagated from the un- 
certainties of FWHM and continuum luminosity mea- 
surements from the spectra; the measurement errors are 
typically ^ 0.3 dex for our sample (see Fig. [U and 
hence are negligible; however, measurement error may 
become important for other samples with low spectral 
quality. 

• statistical error, which is the scatter of virial BH masses 
around RM masses when these virial estimators were 
calibrated against local RM AGN sample; the statis- 
tical error is > 0.3 d ex (e. g., McLure & Jarvis 200 j; 
McLure & DunlopI 12004 IVestergaard & PetersonI 



2006h . which will be taken into account in our 



Bayesian approach. 

• systematic biases, which may result from the virial 
assumption, the usage of RM masses as true masses 
during calibration, the usage of a particular definition 
of line width as the surrogate for the virial velocity, 
the extrapolation of the virial calibrations to high lu- 
minosity/r edhift, as well as other possible systemat- 
ics (e . g., LK i-olik 2001; Collin etal. 2006; Shen et ajJ 
20081: JMarconi et al. 2008; Fine et al. 2008. 2010^ 

I2OO9 



2008; Fine etal. 2008. 
Netzed I2009t iDennev et atl 120091: [Wang et al.1 



Graham et al 1 1201 U iRafiee & Halll l20nbl; ISteinhardtl 

2oTi1k 



We generally neglect systematic biases in the current 
study, as they are poorly understood at present. That 
means we assume on average these virial mass esti- 
mators give unbiased mass estimates (see ^3.2.11 for 
the meaning of "unbiased"). However, we do consider 
a possible luminosity-dep endent bias (e.g., IShen et all 
|2() 08';'Shen & Kellvll2010i) . which we describe in detail 
in 33.2.11 This is not only because luminosity is an ex- 
plicit term in all virial estimators, but also because that 
many studies with virial BH masses are restricted to fi- 
nite luminosity bins or flux-limited samples. Moreover, 
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understanding any potential luminosity-dependent bias 
is crucial to probe the true distribution in the mass- 
luminosity plane. 

This paper is organized as follows. In ^we describe the 
data; we present the traditional binned LF/BHMF in 33. H and 
describe the Bayesian approach in 33.21 We present our model 
results in ^ discuss the results in ^ and conclude in ^ 
Throughout the paper we adopt a flat ACDM cosmology with 
cosmological parameters f^A = 0.7, fio = 0.3, h = 0.7, to match 
most of the recent quasar demographics work. Volume is in 
comoving units unless otherwise stated. We distinguish virial 
masses from true masses with a subscript vir or ^ . Quasar lu- 
minosity is expressed in terms of the rest-frame 2500 A con- 
tinuum luminosity (L = XL\ or / = logL for short), and we 
adopt a constant bolometric correction Cboi = L^oi/L = 5. 

2. THE DATA 

Our parent sample is the SDSS DR7 quasar catalog 
dSchneider et al.l 1201 0'), which contains 105,783 bona fide 
quasars with ;-band absolute magnitude M, < -22 and have 
at least one broad emission line (FWHM > 1000 km s~') or 
have interesting/complex absorption features. Among these 
quasars, about half were tar geted using the fina l quasar tar- 
get algorithm described in lRichards et al.l (12002'). and form a 
homogeneous, statistica l quasar sample (e.g., Richard s et al.l 
l2006tlShen et al.ll2007bh . which we adopt in the current study. 
Quasars in this homogeneous sample are flux-limited to / = 
19.1 below z = 2.9 and to / = 20.2 beyond^. There is also a 
bright limit of / = 15 for SDSS quasar targets, which only be- 
comes important for th e most lum inous quasars at the lowest 
redshift (see Fig. 1 in IShen et al.l l201 1). We have used the 
continuum and emission line /T-corrections in [Richards et al.l 
(12006) to compute the absolute /-band magnitude normal- 
ized at z = 2, Mi[z = 2]. At z < 0.5, host contamination 
beco mes more and rn ore prominent towards lower redshift 
(e.^.. IShen et al.ll201 ll) . so we restrict our sample to z > 0.3. 
Our final sample includes 57,959 quasars at 0.3 < z < 5. 
The sky coverage of this uniform quasar sample is carefully 
determined, using the approach detailed in the appendix in 
IShen et alj (|2007b), to be 6248 degl 

The virial mass estimates and measurement errors for these 

(l20Tll) 



quasars were taken from 



IShen et al.1 (120111) . We refer the 
reader to IShen et all (1201 Ih for details regarding the spectral 
measurements and virial mass estimates. In short, the spectral 
region around each of the three lines (11/3, Mgii, and Civ) is 
fit by a power-law continuum plus iron template"*, and a set 
of Gaussians for the line emission. Narrow line emission is 
modeled for H/3 and Mgii but not for Civ. We use the con- 
tinuum luminosity and line FWHM from the spectral fits to 
compute a virial m ass using one of the fiducial virial cali- 
brations adopted in lShenetaLJ (1201 iL eqns. 5,6,8). > 95% 
of the 57,959 quasars have measurable virial BH masses. 
Fig- [H (right) shows the distribution of measurement errors 
(propagated from the FWHM and continuum luminosity er- 
rors) of these virial mass estimates. The vast majority of 
virial estimates have a measurement error far below 0.3-0.4 
dex, the nominal statistical uncertainty of virial estimators. 

^^ There are a tiny fraction of uniformly-selected quasars targeted by the 
HiZ branch of the target selection algorithm (Richards et al. 2002) at z < 2.9 
dow n to 1 = 20.2. We have rejected these quasars in our flux-limited sample 
(see IShen et alJ2011L for more details). 

* Except for CIV, where we only fit a power-law continuum with no iron 
template applied. 



TABLE 1 

Summary OF zbins 



zbin 


z range 


Nq/N,,, 


M,,Uz = 2] 


H/3 








1 ... 


[0.3,0.5] 


4298/4149 


-22.94 


2 . . . 


[0.5,0.7] 


4206/4027 


-23.84 


Mgll 








3... 


[0.7,0.9] 


3955/3873 


-24.61 


4... 


[0.9,1.1] 


4871/4772 


-25.06 


5... 


[1.1,1.3] 


5872/5789 


-25.39 


6... 


[1.3,1.5] 


5925/5855 


-25.73 


7... 


[1.5,1.7] 


6459/6340 


-25.99 


8... 


[1.7,1.9] 


5839/5566 


-26.29 


CIV 








9... 


[1.9,2.4] 


7761/7545 


-26.83 


10.. 


[2.4,2.9] 


1695/1641 


-27.34 


11.. 


[2.9,3.5] 


4317/4003 


-26.66 


12.. 


[3.5,4.0] 


1830/1666 


-27.00 


13.. 


[4.0,4.5] 


661/518 


-27.36 


14.. 


[4.5,5.0] 


270/152 


-27.45 



Note. — The second column lists the bound- 
aries of each zbin. The third column lists 
the total number of quasars and those with 
measurable virial masses (measurement eiTor < 
0.5 dex) in each zbin. The fourth column 
lists the limiting luminosity in terms of the ab- 
solute /-band magnitude normalized at z = 2 
(Richards et al. 2006), which corresponds to the 
flux limit (i = 19.1 and 20.2 for z < 2.9 and 
z > 2.9) and is estimated at the median redshift 
for each zbin. 



TABLE 2 
Binned DR7 virial BHMF 



logMsHv 

(Mq) 



log'J>(MBH.vii) 



logo-(MBH.vii) 



(Mpc-^ logMaJj^,^) (Mpc-^ logMs[^^■J 



0.4 
0.4 
0.4 



7.50 
7.75 
8.00 



-6.378 
-5.957 
-5.813 



-7.370 
-7.156 
-7.110 



Note. — The full table is available in the electronic version 
of the paper 



Three line estimators were used: H/3 (Vesterg aard & PetersonI 
l2006l z < 0.7): Mgii (Shenetal. 2011. 0.7 < z < 1.9); Civ 
d Vestergaard & PetersonI l2006l ~> 1.9). Virial BH masses 
based on two estimators are smoothly bridged across the di- 
viding redshift, i.e., there is no systematic offset between two 
different estimators. Fig. [T] (left) shows the redshift distribu- 
tion of virial mass estimates in our sample, where the vertical 
dashed lines mark the divisions between two estimators and 
the grid we use to compute the BHMF (see below) is shown 
in gray. We reject objects with a measurement error > 0.5 dex 
in virial mass estimates in computing the BHMF, and we will 
correct for this incompleteness in mass estimates in Sec[3j 

3. THE QUASAR LP AND BHMF 

3.1. The Traditional Approach 

Follow i ng th e common practice in the literature (e.g. . 
Fan et al .' '2001'; 'Richards et al.' '2006'; 'Greene & Hd '2001 
Vesterg aard et al. 2008; Vestergaard & Osmer 200^ 
Schulze & Wisotzki 2010), we use the l/Vmax method 
(e.g.. Schmidt 1968) to estimate the LF and active BHMF: 
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Fig. 1. — Left: Redshift distribution of virial BH masses in our sample. Right. Distribution of measurement eiTors of tlie virial BH mass estimates. The 
vast majority of virial mass estimates have negligible measurement errors compared with the nominal statistical uncertainty of virial BH mass estimators Cvir ~ 
0.3-0.4 dex. 
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the tabulated selection function^ in lRichards et alJ(l2006l) with 
interpolation to estimate Q(L,z), and calculate Vmax for each 
quasar in a redshift-luminosity (virial mass) bin. The binned 
LF, $(M,[z = 2]) = dn/dMi[z = 2] is then: 



$(M,[Z=2]): 



1 



N 



Miiax.j 



-24 -26 -28 -30 



with a Poisson statistical uncertainty 
1 



.DR7 
3DR3 



a($) = 



AMi[z = 2] 



N 

E 



V, 



max, 7 



1/2 



(2) 



(3) 



-24 -26 -28 -30 -24 -26 -28 -30 -24 -26 -28 -30 



Mi [z = 2] 

Fig. 2.— Comparison between the DR3 (Richards et al. 2006") and DR7 
binned LF (this work) for the same luminosity-redshift grid. The DR7 results 
are in good agreement with earlier DR3 results. 



TABLE 3 
Binned DR7LF 



where the summation is over all quasars within a redshift- 
magnitude bin. 

The l/y„ax binned BHMF, ^(Mbumv) = dn/dlogMBHMr, 
is then: 



$(Mb 



r) = 



1 



AlogMBH,v 



N 

■E 



1 



K 



max.j 



M,[z = 2] 



log<I.(M,[.-; = 2]) 
(Mpc^^^mag"') 



log a(M,[z = 2]) 
(Mpc^^^mag"') 



0.4 
0.4 
0.4 



-22.65 
-22.95 
-23.25 



-5.669 
-5.643 
-5.858 



-6.920 
-7.078 
-7.350 



with a Poisson statistical uncertainty 

1 r '^ 
^m = ^rT7ZTi E 

7=1 



AlogMBH., 



V, 



max. 7 



1/2 



(4) 



(5) 



Note. — The full table is available in the electronic 
version of the paper. 



where fl is the sky coverage of our sample, dVc/dz is the dif- 
ferential comoving volume, z^n and z^ax are the minimum 
and maximum redshift within a redshift-luminosity (virial 
mass) bin that is accessible for a quasar with luminosity L, 
and Q{L,z) is the luminosity selection function mapped on 
a two-dimensional grid of luminosity and redshift. We use 



where the summation is over all quasars within a redshift- 
mass bin. 

As a sanity check, we computed the DR7 q uasar luminos- 
ity fun ction (LF) in the same L-z grid as in [Richards et al.l 
(|2006|) . and found it in excellent agreement with the DR3 re- 
sults with smaller statistical error bars (Fig.|2]i. 

To compute the binned LF and BHMF we choose a redshift 
grid (zb ins) that avoids straddling two mass estimators, with 

'' There is no difference in the target selection completeness between the 
uniform DR3 quasars used in Richards et al. (2006) and the uniform DR7 
quasars used here, since the final quasar target algorithm was implemented 
after DRl. 
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boundaries of 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.4, 
2.9, 3.5, 4.0, 4.5, 5.0. Within each of the 14 zbins we use 
a mass grid with a bin size of AlogMeH.vir = 0.25 starting 
from logMBH.vir = 7.375, and a luminosity grid with a bin size 
of AM,[z = 2] = 0.3 stalling from M,[z = 2] = -22.5. Table 
[U summarizes information for each zbin. Fig.[3]shows the 
binned virial BHMF using the l/Vmax technique. Our binned 
virial BHM F results are similar to t he binned virial BHMF 
estimated in lVestergaard et aLl (l2008l) based on DR3 quasars, 
with much better statistics due to increased sample size. 

Two important facts limit the application of the binned 
virial BHMF. First, it is inappropriate to use the selection 
function upon luminosity selection for the BHMF, i.e., BHs 
with instantaneous luminosity fainter than the flux limit of 
the survey will be missed regardless of their masses. As a 
result, the binned BHMF suffers from incompleteness, espe- 
cially at the low-mass end, and the turn-over of the BHMF 
at low masses seen in Fig. [3] is not real. Second, virial BH 
masses are not true masses. Substantial scatter between virial 
mass estimates and the true masses changes the underlying 
BH mass distribut ion, and may lead to signific ant Mal mquist- 
type biases (e.g.. IShenetal.l l2008t iKellv et al. 2009, IJOlOl: 
IShen & Kellyl2010l) . The latter effect is particularly important 
at the high-mass end (where the contamination from intrinsi- 
cally lighter BHs can dominate over the indigenous popula- 
tion) and at high redshift (where the virial BH mass estimator 
switc hes to the more problematic Civ line, e.g., Shen et al. 
|2008|). The Bayesian framework developed in Kellv et ah 
(120091) and described in ^3.2| remedies these issues, and pro- 
vides more reliable estimates for the intrinsic BHMF. 

Bearing in mind the limitations of the binned BHMF, Fig. 
[3] shows a coherent evolution for the most massive (Mbh.vIi }t 
3 X 10^ Mq) BHs: their abundance rises from high redshift 
and reaches maximum around z ^ 2, then decreases towards 
lower redshift. This trend is likely a manifestation of the rise 
and fall of bright quasars seen in the LF, and we will test this 
trend with the Bayesian approach described in Sec |3.2| 

For future comparison purposes only, we tabulated the 
binned virial BHMF in Table|2j but we remind the reader that 
it should be interpreted with caution. We also tabulated the 
binned LF in Table [3] 

3.2. The Bayesian Approach 

As discussed earlier, the causal connection between the LF 
and BHMF naturally requires a determination of the joint dis- 
tribution in the mass-luminosity plane. In doing so, one needs 
to account for selection effects of the flux limit of the sample, 
and to distinguish between virial masses and true masses. The 
best approach is a forward modeling, in which we specify an 
underlying distribution of true masses and luminosities and 
map to the observed mass-luminosity plane by imposing the 
flux limit and relations between virial masses and t rue masses, 
and compare with the observe d distribution (e.g.. lShen et alj 
I2008t IKellv et all l2009l l2010t) . This is a complicated and 
model-dependent problem. Below we first demonstrate our 
best understandings of the relationship between virial masses 
and true masses, then we describe our model parameteriza- 
tions and the implementation of the Bayesian framework. We 
defer the caveats in our model to 



3.2.1. Preliminaries 

Here we describe our modeling of the statistical errors of 
virial mass estimates, under the premise that these FWHM- 
based virial mass estimators on average give the correct mean 



(e.g., see Eqn. 8 below). For clarity, we use p(x\y) to denote 
the conditional probability distribution of quantity x at fixed 
y, andx\y to denote a random value of x at fixed y drawn from 
p(x\y). 

For the local RM AGN sample (which has a dispersion 
of '^ 1 dex in luminosity), single-epoch virial BH mass esti- 
mates were calibrated against RM masses (assumed to be true 
masses) to have the right mean, and a scatter (uncertainty) of 
^ 0.4 dex around RM masses (e.g., McLure & Jarvis 200^ 
iMcLure & Dunlopll2004t IVestergaard & Petersons 2006). To 
account for the effects of the uncertainty in virial BH mass 
estimates, we first must understand the origin of this uncer- 
tainty. It is natural to ascribe this uncertainty to two facts (e.g., 
IShenetalJl2 008: Shen & Kelly 2010): a) luminosity is an im- 
perfect tracer of the BLR size; b) line width is an imperfect 
tracer of the virial velocity. Taken together, some portion of 
the variations in luminosity and line width are independent of 
each other, causing the virial mass estimates to scatter around 
the true value; the remaining portion of the variations in lu- 
minosity and line width cancel with each other, and do not 
contribute to the scatter in the virial mass estimates. 

To better understand this, consider the following example. 
Take a population of A^ BHs with the same true mass m = 
logMfiH, and assuming: a) The FWHM and luminosity follow 
lognormal distributions at this fixed true BH mass; b) a mean 
luminosity -radius (R-L) relation R ex L"^, and a linear mean 
relation between FWHM and the virial velocity v; and c) the 
virial masses are unbiased on average. For this population of 
BHs, the luminosity I = logL of individual object is given by: 



l\m={l)„, + Gi(0\c7',) + Go(OKou) 



(6) 



where l\m is the individual object luminosity at this fixed m, 
G,(/i|cr) is a Gaussian random deviate with mean fi and dis- 
persion a, and (/)„ is the expectation value of luminosity at 
this true BH mass. Similarly we can generate individual line 
width w = log FWHM as: 



w\m= (w}„, + G2(0|cr„,)-0.25Go(0|crcoiT) 



(7) 



where {w)m is the expectation value of line width at this true 
BH mass. The individual virial mass estimate me = logMBH.vir 
at this fixed Mbh is then 



me\m = m + 0.5Gi(0\a'i) + 2G2(0\a„) 



(8) 



which implies that the virial BH mass estimates follow a 
lognormal distribution around the correct mean (i.e., m), 
but have a lognormal scatter (virial uncertainty) a^ir = 
y/(0.5a')^ + (2a„)^ around the mean (e.g.. IShen etall 120081; 
IShen & Kellyll2010l) . The Go terms of variation in luminosity 
and FWHM exactly cancel with each other and do not con- 
tribute to the virial uncertainty, and were referred to as the 
"correlated variations" in FWHM and luminosity in the above 
papers; while the Gi and G2 terms were referred to as the 
"uncorrected variations" in FWHM and luminosity, and they 
contribut e to the virial uncertainty in quadratic sum. The ap- 
proach in lKelly et all ('2009', '20 10*) implicitly assumed a', = 
while the approach in Shen et al. (2008) and Shen & Kellvl 
(I2OIOI) is to set (Tcon = and consider non-zero a'l- The latter 
choice is motivated by the fact that the observed distribution 
of FWHM for SDSS quasar samples is already narrow (dis- 
persion < 0. 15 dex) and the premise that the virial uncertainty 
(Tvii should be no less than ^ 0.3 dex. 

Physically cr'i is unlikely to be zero. If this were true, it 
would imply that single-epoch luminosity is an unbiased in- 
dicator for the instantaneous BLR radius at fixed BH mass. 
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TABLE 4 
Model LP. BHMF and Eddington ratio function 
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logZ. 
(ergs-') 


(Mpc--'dcx- 
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logMBH 
(Mq) 


*0 *+ 

{Mpt--'(lcx-' 


*- 


*0 *+ , 
(Mpc-'dcx-' 


*- 


log A 


«o 


{Mpt-3jcx-' 


) 


(Mpc^^dcx"' 


<^- 


0.4 
0.4 
0.4 


-16.904 
-16.979 
-17.054 


42.00 
42.03 
42.06 


-6.299 -6.002 
-6.240 -5.950 
-6.182 -5.900 


-6.639 
-6.572 
-6.506 


6.000 
6.025 
6.050 


-10.684 -8.179 
-10.550 -8.114 
-10.415 -8.049 


-12.383 
-12.201 
-12.023 


-19.804 -17.507 
-19.521 -17.287 
-19.240 -17.068 


-21.693 
-21.363 
-21.038 


-1.000 
-3.975 
-3.950 


-9.482 
-9.368 
-9.256 


-9.108 
-9.002 
-8.896 


-10.225 
-10.097 
-9.970 


-17.622 -16.349 
-17.397 -16.152 
-17.175 -15.955 


-19.266 
-18.999 
-18.736 



Note. — The full table is available in the electronic version of the 



While in practice it is more natural to expect there are uncor- 
related random scatter in both L and R, indicating a stochastic 
term in addition to the deterministic term (when predicting R 
with L), which will lead to biased estimates for R at fixed L. 
The sources of this stochastic term may include: a) the con- 
tinuum luminosity variation and response of the BLR is not 
synchronized; b) individual quasars have different BLR prop- 
erties; c) optical-UV continuum luminosity is not as tightly 
connected to the BLR as the ionizing luminosity. Further- 
more, even if single-epoch luminosity were an unbiased in- 
dicator of the instantaneous BLR radius, certain line width 



indicators (such as FWHM) might still not response to all the 
variations in luminosity. For example, consider a single BH 
where its luminosity varies (and its BLR radius varies instan- 
taneously following a perfect R-L relation), and suppose that 
the broad line is composed of a non-virialized component and 
a virialized component. When luminosity increases (BLR ex- 
pands), the virialized component reduces line width, but the 
non-virialized component may increase line width if it is orig- 
inated from a radiatively driven wind (in the case of Civ, more 
blueshifted Civ tends to have a larger FWHM, e.g., Shen et al. 
2008): the combined line FWHM may not change due to the 
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two opposite effects. Therefore in this case although luminos- 
ity is tracing the BLR size perfectly, some variation in lumi- 
nosity is not compensated by variations in FWHM and should 
be counted as the uncorrected variation a\. 

A non-zero a\ implies that the distribution of virial mass 
estimates at fixed true mass and fixed luminosity, p(me\m,l), 
is different from the distribution of virial mass estimates at 
fixed true mass, p{me\m). In the extreme case where Ccorr = 
0, i.e., FWHM does not change in response to variations in 
luminosity at all, we have 

m,|m,/ = m + O.5(/-(0m) + 2G2(Olcr„) . (9) 

Hence not only is the distribution p(me\m,l) narrower than 
pinie \m), but also the expectation value of nig is biased from 
the true BH mass for any fixed luminosities except for / = (/),„. 
Now consider a more general form of the luminosity distri- 
bution at fixed true mass and the actual slope in the observed 
luminosity-radius relation, we can parameterize the distribu- 
tion of p(me\m, I) as: 



\m,l ■■ 



i + l3(l-{l)„) + e„i , 



(10) 



where again (Z)„, is the expectation value of luminosity at 
fixed true mass, e,„/ is a random deviate with zero mean and 
dispersion cr„,/, denoting the scatter of virial mass estimates 
affixed true mass and fixed luminosity, and the error slope /3 
describes the level of luminosity-dependent mass bias at fixed 
true mass and luminosity. Both (3 and e„,/ are to be constrained 
by our data. Eqn. (fTOl i implies that the variance of mass esti- 
mates at fixed true mass and luminosity is reduced to: 



Var(OTe|OT,0 = Var(OT(.|OT)(l -p ) 



(11) 



where p^ = /3^Var(/|m)/Var(me|m), and Var(...) refers to the 
variance of a distribution. The formal uncertainty of the virial 
mass estimator is then 



(Tvir = ■\/Var(m7|ffj)= y Vai(me\m, I) + /S^Waiillm) . (12) 

If we assume a single log-normal distribution for p(l\m) and 
emi (with a dispersion ami), the above equation reduces to 



(Tvir = 



ml 



+/3V/ 



2^2 



(13) 



Note that here cr/ is the total dispersion in logL at fixed mass, 
rather than the portion cr,' that is not responded by FWHM as 
in Eqns. (|6]l and (|8]l. 

Eqn. ( [Tol l is a rather generic form that describes the rela- 
tion between virial masses and true masses and the possible 
luminosity-dependent bias in virial masses^, and is one of the 
basic equations in our Bayesian approach. The value of (3 de- 
pends on the relative contributions from cr,' and cTcon in the lu- 
minosity dispersion at fixed mass. Under the assumption that 
the mean R-L relation and a linear mean relation between 
FWHM and virial velocity are correct as in the adopted virial 
estimators, a non-zero cr,' leads to a positive /3. If the value 
of /3 approaches the slope in the adopted mean R-L relation, 
then it suggests either luminosity or FWHM is a poor indica- 
tor for BLR size or virial velocity over the narrow dynamical 
range at fixed true BH mass (although they could still be rea- 
sonable indicators for large dynamical ranges in mass and lu- 
minosity). On the other hand, if /3 is small, then it means lumi- 
nosity and FWHM vary in concordance even over the narrow 

* One can work out a similai' equation for the distribution of m^ at fixed m 
and FWHM w, p{ine\m,w) = in + 0' (w — {w) „,) + ^mw ■ A non-zero cr,,. in Eqn. 
(7) will lead to a non-zero /3' and p{me\m, w) ^ p(me \m). However, this is of 
little practical value since virial masses are never binned in FWHM. 



dynamical range at fixed true BH mass, and hence are good 
indicators for BLR radius and virial velocity. /3 = represents 
the extreme situation where FWHM responds to all the varia- 
tion in luminosity at fixed true mass (plus additional scatter in 
FWHM), and no bias in virial masses is incurred when lumi- 
nosity deviates from (Z),„. /3 = is generally assumed in most 
studies with virial BH masses. 

We also note that one advantage of using Eqn. ( fTOl i is that 
it does not rely on the assumption that the mean R-L rela- 
tion and the linear mean relation between FWHM and virial 
velocity used in these virial estimators are correct. In other 
words, if the virial estimators adopted in this work used incor- 
rect forms for the mean R-L relation and the mean relation 
between FWHM and virial velocity, then a negative value of 
(3 may be needed to correct nig at fixed m and /. Of particular 
interest here is whether or not ra diation pressure is irn portant 
in the dynamics of the BLR (e.g.. lMarconi et an2008l) . which 
will indicate a negative (3 based on the virial masses that have 
not been corrected for radiation pressure. We will test if a 
negative (3 is required to model the observed distribution in 
our Bayesian approach. 

Is there any indication for a non-zero (3 from the reverber- 
ation mapping AGN sample ? There are only ^ 3 dozens of 
RM AGNs and we do not have enough objects with the same 
BH mass to test source-by-source variations. Nevertheless, 
we can still test the luminosity-dependent bias using repeated 
spectra for the same object when its luminosity changes a sig- 
nificant amount over time. NGC 5548 is the most frequently 
monitored RM AGN (H/3 only), and has been observed in 
different luminosity states with a spread of ^ 0.5 dex in lu- 
minosity (e.g., Peterson et al. 2004; Bentz et al. 2009b), thus 
provides an ideal test case for single-source variations. 

In Fig. m we show the H/3 virial product for NGC 5548, 
computed using the continuum luminosity and line width 
measured at different luminosity states in each monitoring 
period, as a function of continuum luminosity . The spec- 
tral measurements were taken from ICoUin et al.l 12006), and 
we have corrected the continuum l uminosity for host st arlight 
using the correction provided by iBentz et al] (l2009ah . The 
line widths were measured from both the mean and rms spec- 
tra^ for each monitoring period. The left and right panels 
of Fig. |4] show the virial product computed using FWHM 
and line dispersion, respectively, and its scaling with lumi- 
no sity is the same as in the viria l mass estimators provided 
bv IVestergaard & Petersoiil (|2006|) . The FWHM-based virial 
product shows an average trend of increasing with luminos- 
ity, which means that FWHM does not fully response to 
the variations in luminosity, leading to a positive bias in the 
virial product (and thus in the virial mass estimate). This 
trend seems to be shghtly weaker when using line disper- 
sion instead. A lin ear reg ression analysis using the Bayesian 
method of Heliv] ^W) yields: /3 - 0.65 ± 0.27 (FWHM, 
mean); ;3 ~ 0.51 ±0.34 (FWHM, rms); /3 - 0.20 ±0.30 {ar,ne, 
mean); (3 ^ 0.45 ± 0.29 (criine, rms), where uncertainties are 
Icr. While it is inconclusive based on this single object, there 
is some indication that a positive j3 is favored, especially for 
the virial product based on FWHM from the mean spectra, 
which is the closest to that based on FWHM from single- 
epoch spectra. It would be important to test this for more 

' Strictly speaking, for single-epoch virial mass estimates, neither the 
mean nor rms spectra are available. However, the spectral variability dur- 
ing each monitoring period is small enough such that the mean spectrum is 
close to single-epoch spectra within this period. 
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Fig. 4. — The dependenc e of the virial prod uct computed from luminosity and line width as a function of luminosity, for a single object NGC 5548 and for 
H/3 only. The data are froml CoUin et al.l J2006I) . and are based on both mean and rms spectra during each monitoring period. Error bars represent measurement 
errors. The en'or bars in luminosity have been omitted in the plot for clarity. The continuum luminosity has been corrected for host starlight using the coiTection 
provided by Bentz et alJ 1 2009a). the black and blue dashed lines are the best linear-regression fits using the Bayesian method of Kellv 1 2007), for measurements 
based on mean and rms spectra, respectively. Left: virial product based on FWHM; the data point for Year 5 (JD 48954-49255) based on the rms spectrum has 
been suppressed due to problematic measurements (e.g.. , Peterson et al. 2004; Collin et al. 2006). Right: virial product based on line dispersion CTHne- 
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Fig. 5. — Model parameters for zbin2. Shown here are the posterior dis- 
tributions of some model parameters and derived quantities. From top-left in 
clockwise order: the dispersion in mass estimates at fixed true mass and lu- 
minosity, cr^i; the error slope /3; the slope in the mean (true) mass-luminosity 
relation for our Eddington ratio model, aj ; the dispersion in Eddington ratios 
for all broad-line quasars, CT(log A) where A = Lbol/^Edd^ 'he mean Edding- 
ton ratio for all broad-line quasars, (log A); the dispersion in luminosity at 
fixed true mass in our Eddington ratio model, ct; . 

objects with repeated spectra, and for Mgii and Civ as well. 

To summarize, because luminosity is an explicit term in 
virial mass estimators, these virial mass estimates are no 
longer independent (and unbiased) estimates of true masses 
when restricted to a narrow luminosity range or a flux-limited 
sample, for cases where fi ^ Q. Our view of the uncer- 
tainties in these virial mass estimates (i.e., the scatter in 
p(me\m), as determined in the calib rations against t h e RM 
AG Ns) is thus different from that in iKollmeier et all (l2006h 
and lSteinhardt & Elvisl (l2010bh . 

3.2.2. Implementing the Bayesian Framework 



45.0 45.5 46.0 
log L [erg s"'l 



7.0 7.5 8.0 8.5 9.0 9.5 lO.O 10.5 
log M„H.„ [MJ 



Fig. 6. — Posterior checks for zbin2. The solid black histograms shows 
the observed distributions. The red points and error bars are median results 
and uncertainties from our simulated samples using 500 random draws from 
the posterior distributions. The top and bottom panels show the histograms 
in linear and logarithmic scales, respectively. 



Now we proceed to describe our model setup and the im- 
plementation of the Bayesian framework. Below we describe 
the basics of our model approach. M ore details regardin g the 
Bayesian approac h can be found in iKellv et all (|2009() and 
lKellvetalJ(l2OT0l) . 

a. The BHMF and lu minosity distribution model. As in 
iKelly et alJ (l2009l IToiOi) . we use a mixture of log- 
normal distributions as our model for the true BHMF, 
and a single log-normal luminosity (Eddington ratio) 
distribution at fixed true BH mass. The mixture of log- 
normals is flexible enough to capture the basic shape 
of any physical BHMF, and greatly simplifies the com- 
putation as many integrations can be done analytically. 
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The model true BHMF reads 

1 K 



ci>M=A^(^) y: 

k=\ 



dz 



TTt 



iTial 



:exp 



(m-^k) 



2al 



, (14) 



where m = logMfiH, N is the total number of quasars, 
/Xjt and ai! are the mean and dispersion of the A;th Gaus- 
sian, and J2k=i "'* = 1- We use K = 3 log-normals to 
describe the BHMF, as we do not find significant differ- 
ence when increasing the number of log-normals used. 
The luminosity distribution at fixed BH mass is mod- 
eled as 



p(l\m) ■■ 



1 



Inaj 



:exp 



[l-aQ-ai(m-9)y 

2^f 



(15) 



where / = logL, ao and ai describe a mass-dependent 
mean luminosity, and ai is the scatter in luminosity at 
fixed mass. The LP is therefore 



(16) 



$(0= / <^{m)p(l\m)dm . 



The virial mass prescription. We assume that virial 
masses are unbiased whe n aver aged over luminosity at 
fixed true mass (i.e., Eqns. 18110b . and we generate virial 
masses at fixed true mass and luminosity according to 
Eqn. ( fTOl ), assuming a single Gaussian (with dispersion 
(T„,/) to describe the scatter e,„/ at fixed mass and lumi- 
nosity. 

The redshift distribution. Because the redshift bins are 
narrow, we approximate the distribution of redshifts 
across the bin as a power-law, where the power-law in- 
dex 7 is a free parameter: 



p(z\j) ■■ 



(1+7)^'' 



1+7 
^max " 



1+7 ' 



(17) 



Here, Zm^x and Zmin define the upper and lower boundary 
of the redshift bin, respectively. 

d. The posterior distribution p(9\mg, I, z) is 

N 

p{9\me,l,z)^p{e)[p{I=im-''Y[p(m,^i,li,Zi\e) (18) 

where 0(7rA.,/i/.,crA.,ao,ai,cr/,/3,iT,„/,7) is the set of 
model parameters, A^ is the total number of quasars, 
p{d) is the prior on 6, p(I = l\9) is the probability as a 
function of 9 that a broad-line quasar is included in the 
flux-limited SDSS quasar sample, and the likelihood 
function p(mp,-,/,,z,|6') is determined by Eqns. (fTOl l. 
dUli, dBJ and dnii. 

In this work, we derive the continuum luminosity at 2500 A, 
I = logL, from the /-band magnitude according to the pre- 
scription given in Richards et al. (2006). This is a departure 
from the approach taken by Kellv et al. (201 0|), who used th e 
continuum luminosity esti mated bv.Vestergaard et aP (l2008l) . 



However, as explained in iKelly et aP (1201 Ol) . because the 



SDSS selection function is in terms of the /-band magnitude, 
they had to assume a model for the distribution of / at fixed 
luminosity, and then calculate the selection function in terms 
of luminosity by averaging over this model distribution. They 
note that this approach can lead to instability in the estimated 
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Fig. 7. — Posterior checks for zbin2 in the mass-luminosity plane above 
the flux limit. The two black lines indicate Eddington ratios of 0.01 and 1. 
The black and red contours are for the observed and simulated distributions 
using virial masses, and the blue dashed contour shows the simulated distri- 
bution with true masses. Our model fits the observed distribution well, and 
the distribution using true masses is different from that using virial masses. 

mass function in mass bins that are severely incomplete, as 
small errors in the selection function can lead to large de- 
viations in pil = l\6), which appears in the denominator in 
Equation (fTsT l. Instead, we use the luminosity derived from 
the /-band magnitude to ameliorate this effect, as the selec- 
tion function is calculated in terms of /. 

It is necessary to impose some prior constraints on /3 and 
(T,„i based on the reverberation mapping data set, as these pa- 
rameters are degenerate with some of the other parameters. 
Unlike most previous work, we do not fix the values of (3 and 
cr,„/ to, say, /? = and cr,„/ = 0.4 dex, but use a prior distribution 
which incorporates our uncertainty in these parameters. This 
uncertainty will be reflected in the probability distribution of 
the mass function, given the SDSS DR7 data se t . We a pplied 
the Bayesian linear regression method of iKellv I (l2007h to the 
reverberation mapping sample in order to estimate the prob- 
ability distr ibution of /? and a,ni based on this sample. The 
method of Kellv I ( l2007h incorporates the measurement errors 
in the mass estimates, which is important when estimating 
the amplitude of the scatter in the mass estimates. We set 
(/)„, equal to the mean luminosity for the reverberation map- 
ping sample^. We used the values o f RM black ho l e mas s 
(assumed to be true masses) given bv ' Peterson et al.l (|2004|) . 
For the H/3 calibration, we used the value of 5 1 00 A luminos- 
ity givenjn Bentz et al. (2009a) and value of FWHM given 
in IVesterpaard & Petersoni 62006) ; when there were multiple 
measurements, we averaged them together. For Ci v we used 
the values given in IVestergaard & Petersoni (2006). For H/? 
we found that (3 = 0.16±0.1, and that the posterior distribu- 
tion for (tI^i is well described by a scaled inverse x^ distribu- 
tion with v Ri2Q degrees of freedom and scale parameter 5^ = 
0. 1. For Civ we found that /3 = 0. 15 ±0. 14, i^« 20, 5^ = 0. 17. 
For Civ this is similar to the usually quoted scatter in the mass 
estimates of s « 0.4 dex, but the scatter is less for H/3, s « 0.3 
dex. This reduced uncertainty is likely because we have used 

** Ideally one should use a population of obj ects with the same true BH 
mass to perform linear regression using Eqn. (To), which is not available 
given the limited size of the RM sample. So instead we use the whole RM 
sample. Nevertheless, the derived prior constraint on /3 is weak and our re- 
sults are insensitive on the prior on /3. 



10 



SHEN & KELLY 



45.5 



44.0 



M 43.5 



43.0 



42.5 




Fig. 8. — The simulated mass-luminosity plane for zbin2, which extends 
below the flux limit (the black horizontal line). The red contour is the dis- 
tribution based on true BH masses, and is determined by our model BHMF 
and Eddington ratio model. The black contour is the distribution based on 
virial BH masses. The flux limit only selects the most luminous object into 
our sample, and the distribution based on virial BH masses is flatter than the 
one based on true masses due both to the scatter cr,,,/ and a non-zero /3 (see 
Eqn.fTol. 

5100 A l uminosity values wh ich are corrected for host galaxy 
starUght (iBentz et al.ll2009al) . 

Based on the reverberation mapping results, we impose a 
Cauchy prior distribution on (3 with mean and scale param- 
eters equal to those derived from the reverberation mapping 
data, and we impose a scaled-inverse x^ prior distribution on 
(j-^i with 1/ = 15 degrees of freedom and scale parameter set to 
that from the reverberation mapping sample. For Mgii we use 
the values derived for H/3 since the single-epoch virial masses 
base d on the two hnes seem to correlate with each other well 
(e.g.J Shen et al.l201 Ih . We have used a Cauchy prior because 
it has significantly more probability in the tails than the usual 
Gaussian distribution, making our prior assumptions more ro- 
bust. Similarly, we reduced the degrees of freedom for the 
prior on ct,^^, compared to the reverberation mapping sample 
in order to increase our prior uncertainty on ami ■ This choice 
of prior assumes an uncerta inty on ami of ^ 20%. 

As in lKellvetall (120091) and iKellv et aTl (120101) . we use a 
Markov Chain Monte Carlo (MCMC) sampler algorithm to 
obtain random draws of 9 according to Eqn. ( fTsT i and thus the 
posterior distribution of model parameters given the observed 
data in the virial mass-luminosity plane. Our MCMC sam- 
pler employs a combination of Metropolis-Hastings updates 
with p arall el tem pering. The reader is referred to Kel ly et alj 
(|2009) and Kelly et al. (2010) for further details. 

Different from .Kelly et ak C2009.) and iKellv etakl (1201 Ol) . 
we model the data in individual redshift bins instead of for the 
whole sample. We treat each redshift bin as an independent 
data set. This allows us to explore possible redshift evolution 
of our model parameters, as well as the sensitivity on changes 
in the detection luminosity threshold and the specific virial 
mass estimator used. The caveat is that the constraints are 
generally weaker given less data points in each bin, and that 
the constrained parameters do not necessarily vary smoothly 
across adjacent redshift bins. 

4. RESULTS OF THE BAYESIAN APPROACH 

4. 1 . zbin2 as an example 
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Fig. 9. — Model LF (top panel) and BHMF (bottom panel) for zbin2. 
Th e data points and en'or bars are the binned LF and virial BHMF estimated 
in jl3.ll The color shaded regions are the 68% percentile range from our 
model LF and BHMF. In the bottom panel, the green shaded region is for the 
detectable (i.e., above the flux limit) true BHMF, and the magenta one is for 
all the broad-line quasars. The turnover of the magenta line below ~ 10* Mq 
is a feature constrained by the data and our model, i.e., if there were more 
lower mass BHs, it would be difficult to fit the observed distribution in the 
mass-luminosity plane (cf. Fig. [8). 



We use zbin2 as an example to demonstrate the infor- 
mation that we can retrieve from the posterior distributions. 
This bin uses the most reliable H/3 line to estimate virial BH 
masses; it also has negligible host galaxy contamination com- 
pared with zbinl. Therefore the constraints for this bin are 
expected to be the most robust. 

Fig.|5]shows the posterior distributions of our model param- 
eters for zbin2. These parameters are tightly constrained, 
although degeneracy does exist among these parameters. 

Fig. |6] presents the posterior checks of our model against 
the data. The black histograms show the distribution of ob- 
served luminosities and virial masses. The points and error 
bars are the median and 68% percentile for simulated sam- 
ples generated using random draws from the posterior distri- 
bution. This ensures that our model reproduces the observed 
luminosity and virial mass distributions. Fig. |7] further shows 
the comparison between model prediction and data in the two- 
dimensional mass-luminosity plane, where the model is the 
one that has the maximum posterior probability. The black 
and red contours show the observed and model-predicted 
joint-distributions of luminosity and virial mass in our sam- 
ple, while the blue contours show that for the true masses. 
The true masses are scattered and biased according to Eqn. 
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Fig. 10. — Eddington ratio function for zbin2. As in the bottom panel of 
Fig.|9] the green and magenta shaded regions are the 68% percentile range for 
the detectable and all broad-line quasars, calculated using our model posterior 
distributions. The top and bottom panels show the Eddington ratio function 
in logarithmic and lineal' scales, respectively. While most broad-Hne quasars 
have a mean Eddington ratio of ~ 0.02 for this redshift, the detected ones 
have a higher mean Eddington ratio of ~ 0. 1 due to the Malmquist-type bias. 

( fTOb . For this particular bin the luminosity-dependent bias be- 
tween virial and true masses is only moderate (/? ^ 0.15 as 
seen in Fig.|5]). 

Fig. [8] summarizes the relation between luminosity and BH 
mass, the effects of the flux limit, and the difference between 
virial masses and true masses in the mass-luminosity plane. 
Black and red contours indicate the simulated distributions 
using virial masses and true masses, respectively. Quasars 
in this zbin peak around 10^ Mq with sub-Eddington ratios. 
This turn-over at low mass is a feature as constrained by our 
model. If the low-mass end BHMF continued to rise, there 
would be more smaller BHs with high Eddington ratio being 
scattered into our sample, and it would be difficult to fit the 
observed distribution above the detection threshold (the black 
horizontal line). However, since the low-mass end is only con- 
strained by a few objects close to the flux limit, the location of 
the turn-over will most likely depend on our model assump- 
tions, such as the mixed-Gaussian model for the underlying 
BHMF and the simple log-normal Eddington ratio distribu- 
tion at fixed true mass; the assumption that the error distri- 
bution of virial masses is Gaussian may also affect the result. 
To fully tackle these issues, deeper data is needed to provide 
more stringent constraints at the low-mass end, and different 
models (such as a mixture of power-laws for the true BHMF), 
albeit more computationally challenging, are required to test 



the robustness of this turn-over Therefore we do not claim a 
robust detection of a turn-over in the true BHMF in this work, 
but simply note the possibility of such a turn-over, and point 
out that the turn-over at larger masses seen in the flux- limited 
BHMF is a selection effect (see below). 

The luminosity at fixed BH mass has a substantial scat- 
ter. However, the SDSS flux limit only picks up the lu- 
minous objects. Since there are more low-mass BHs with 
high Eddington ratios being scattered into our sample, the 
Eddington ratio distribution i s subject to a Malmqui st-type 
bias (e.g., |E ddington| |1913|; iMalmquistI [1922; LaueretalJ 
I2007tisiiene t al. 2008; Kefly et al. 2009, 2010). Furthermore, 
these BH masses above the detection threshold are estimated 
as virial BH masses according to Eqn. ( fTOl l. which further 
stretches the distribution horizontally in the mass-luminosity 
plane. The flux limit, the scatter in m^ at fixed m (cr,,,/), and 
the luminosity-dependent bias (inferred by the value of /?), 
have changed the distribution in the mass-luminosity plane 
substantially. There is much w eaker evidence for a "sub- 
Eddington boundary" claimed by'S teinhardt & Elvisl (1201 Oah 
if we use true masses instead of virial masses, and there is no 
need to inv oke alternative virial calibrations to remove this 
boundary (.Rafiee & Halll 1201 lab . However, such a bound- 
ary may exist if the current versions of virial estimators are 
systematically biased, or if our model parameterization is in- 
appropriate. A more detailed investigation using a different 
model parameterization of the mass-luminosity plane will be 
presented elsewhere (Kelly & Shen, in preparation; see dis- 
cussions in 35.3l l. 

Fig.|9]shows the model LF (upper) and BHMF (bottom) and 
comparisons with data. The shaded curves indicate the 68% 
percentile range. The LF is well constrained even when ex- 
trapolated to fainter luminosities than probed by the SDSS 
sample. The true BHMF is plotted in magenta. It turns 
over below ^ IQ^ Mq, as already noted in Fig. [8] The green 
shaded curve indicate the detected BHMF, i.e., the population 
of quasars that have instantaneous luminosities above the flux 
limit of the SDSS sample. The flux limit causes significant se- 
lection incompleteness in terms of BH mass, which becomes 
worse towards smaller BHs. As a consequence, the turn-over 
in the detected BHMF shifts to a larger BH mass. We also 
over-plotted the binned virial BHMF for our flux-limited sam- 
ple in open circles. The Poisson errors in the virial BHMF 
are not meaningful, and significantly underestimate the un- 
certainty in the BHMF. The shape of the virial BHMF also 
differs from the green curve, due to the difference between 
true and virial BH masses''. 

Fig.[TO]shows the Eddington ratio function of all broad line 
quasars in magenta, and of those detectable quasars in green. 
Most quasars are accreting below Eddington, and the Edding- 
ton ratio distribution for the detected quasars is biased high 
due to selection effect. There is a significant dispersion (~ 0.5 
dex) in Eddington ratios at fixed BH mass and for the entire 
quasar population. 

4.2. Results in all zbins 
Fig, nn shows the 68% percentile range of parameters a,„i, 
13, ai, Qfo, Oil, and Uvir = \/<^mi + (^^f^f- The green, cyan and 
red colors indicate the H/3, Mgii and Civ mass estimators, re- 
spectively. The constraints on some parameters (such as /3, ai 

' We note that in general dm,, j^dm, so the integrated area under the green 
curve and the virial BHMF does not necessarily agree. But the total number 
of "detected" quasars in our model is the same as observed. 
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Fig. 1 1. — Model parameters for all 14 zbins. The green, cyan and red colons are for the H/3, Mgll and CIV respectively. The vertical extent of each segment 
indicate the 68% percentile range from our posterior distributions. 



and cki) are poorer than the others, and there is generally de- 
generacy among these parameters. Also, due to the changes of 
detection luminosity threshold with redshift and the complex- 
ity of the multi-dimensional parameter space, adjacent red- 
shift bins do not necessarily make smooth transition, even if 
we stick to one mass estimator 

One purpose of our Bayesian approach is to explore if 
there is evidence for a luminosity-dependent bias in virial 
BH masses through the parameter (3 in Eqn. (fTol i. to be con- 
strained by the data. We found that for H/3 there is generally 
no need for such a bias. However, for Mgii and Civ, there is 
some evidence that the best models favor a positive /3 value 
between 0.2 and 0.6. Because of this positive /3 value, the av- 
erage virial BH masses are biased high by a factor of a few at 
z > 0.7. Including this luminosity-dependent bias, the uncer- 
tainty in virial mass estimates is a^n- > 0.3 dex and larger than 
the scatter cr,,,/ at fixed luminosity and true mass. We note that 
this positive luminosity bias is in t he opposite sense to the bias 
proposed by lMarconi et al.l (l2008l) . who argue that these virial 
estimators do not include the effect of radiation pressure and 
hence are biased low at high luminosities. This would argue 
for a negative value of (3, which is not seen in our MCMC 
results. 

The dispersion in luminosity at fixed BH mass is con- 
strained to be 0.2 ^ (Ti < 0.6 with a median value of ^ 0.4 
dex. The slope in the mean luminosity and mass relation 
is constrained to be 0.6 < ai < 1-2, and there is evidence 
that the normalization ao increases with redshift. These co n- 
straints are weaker than earlier work (iKelly et al.l2 009. 2010), 
due to the additional freedom introduced by the luminosity- 
dependent bias and the fact that we are now modeling the data 
in individual redshift bins independently. 

Finally, as a sanity check. Fig. [12] shows the compari- 
son between model and data for all 14 zbins in the mass- 



luminosity plane, where again the model is the one that has the 
maximum posterior probability. In all zbins the observed 
distribution is reproduced. The true mass distribution is gen- 
erally different, but the level of this difference varies from bin 
to bin. The latter reflects the large uncertainty of determining 
the relation between true mass and virial mass (i.e., Eqn.fTOli 
when the luminosity threshold and line estimator change with 
redshift. 

4.3. The LF and the BHMF 

Given the full posterior distributions of model parameters, 
we now proceed to compute the model LF and BHMF using 
Eqn. (fl4]l-(fT6]l. Instead of using the median and percentile 
of model parameters, we estimate the median and 68% per- 
centile of the LF and BHMF from individual LF and BHMF 
realizations. We tabulate the model LF and BHMF in Table|4] 

Fig.[T3]compares our model LF with the observed LF. The 
black data points are the binned DR7 LF, and other color- 
coded data points are from various deeper quasar surveys 
(Bongiorno et al. '20071 ICroom et al.ll2009t llkeda et alllloilt 
Glikman et al. 2011). Since the LF data for deeper surveys 
were not necessarily measured in exactly the same redshift bin 
as our grid, we plot them in the closest bins possible and we 
indicate their redshifts in the corresponding colors. The green 
lines are our model LF and the 68% percentile range. The 
conversions betw een different niagnitudes and M,[z = 2] are 
given below (e.g.. [Richards et al.ll2006t ICroom et"al]|2009l) : 

/1450A\ 
M,[z = 2]=Mi450-2.5a^log — -^ -2.5(1 + a^)log3 

V7471A/ 

=MB(vega)- 0.804 

/4670A\ 



--Mg[z = 2]-2.5a^\og 



V7471A/ 
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Fig. 12. — Posterior checks for all zbins in the mass-luminosity plane above the detection threshold. The black contours are for the observed quasars with 
virial masses; the red and blue contours are for the simulated quasars using the model realization with the maximum posterior probability, with virial masses 
and true masses, respectively. The two straight lines indicate Eddington ratios A = 0.01 and 1. The redshift is marked in each panel. Our model reproduces the 
observed distributions well, and demonstrates the difference between true masses and virial masses. 

eters do, such as the mean and dispersion of Eddington ratios 
at fixed true mass. There are, however, some notable discrep- 
ancies in the highest redshift bins between our model extrapo- 
lation and the observed faint-end LF. In particular, our model 
extrapolation is unable to match the faint-end LF results in 
Glikman et al. (2011) and Ikeda et al. (201 1) at z - 4.25. This 
is most likely caused by the failure of our model extrapola- 
tion due to the much poorer statistics and systematics with 
the Civ-based virial masses in the two highest redshift bins. 

Fig.[T4lshows our model BHMF in broad-line quasars. The 
magenta shaded region indicates the 68% percentile of the 
true BHMF, and the green shaded region indicates the 68% 
percentile of the portion of BHMF that can be detected in the 
flux-limited SDSS sample. The data points are the binned 
virial BHMF measured in EU The flux Hmit of SDSS 
greatly reduces the abundance of active BHs towards the low- 
mass end. At the same time, the shape and peak BH mass 



where d= 10 pc, c is the speed of light, A = 2500 A and we 
assume a continuum slope a^ = -0.5. 

Our model LF is constrained by the SDSS quasars alone, 
but it also provides reasonable prediction when extrapolated 
to ^ 3 magnitudes fainter. This means our BHMF and Ed- 
dington ratio model is reasonable when extrapolating not too 
far be yond the regime probed by SDSS quasars. ISchafed 
(l2007h used a purely statistical technique to extrapolate the 
[Richards et all (l2006h DR3 LF to fainter magnitudes. Com- 
pared with his method, our extrapolation is more physically 
grounded: it is constrained by a physical model describing the 
underlying BHMF and Eddington ratio distribution; although 
the mixed-Gaussiaon parametrization for the true BHMF does 
not have any particular physical meaning, some model param- 



14 



SHEN & KELLY 



cn 
D 



I 



O 
Q_ 



CN 

II 

N 

^" 




22-24-26-28-30 -22-24-26-28-30 

■SDSS-DR7 D2SLAQ-SDSS 
oVVDS 
aCOSMOS 
aNDWFS-DLS 



-22-24-26-28-30 -22-24-26-28-30 



M, [z = 2] 



Fig. 13. — M odel LF and comparisons with data. The points represent LF measurements from different surveys: SDSS DR7 (black, this work); 2SLAQ-SDSS 
fblue. lCroonT et al. 2009); VVDS (magenta, Bongiorno et al. 2007); COSMOS (cyan, Ikeda et al. 2011); NDWFS-DLS (red, Glikman et al. 2011). The green 
solid and dashed Hnes indicate the median and 68% percentile of our model LF. The constraint in a zbin is generally better if there are more quasars in that bin. 
The redshift of each zbin is marked in the upper-right corner, and the redshifts for other LF data are marked in the lower- left comer in the corresponding colors. 



are generally different for the true BHMF and for the ob- 
served virial mass function, which is caused by the difference 
between true masses and virial masses (see Eqn. [TOb . The 
model BHMF (for all active BHs or detected BHs) has a much 
larger uncertainty than indicated by the Poisson errors associ- 
ated with the binned virial mass function. This large uncer- 
tainty is mostly caused by the flexibility of our model and the 
poorly constrained luminosity-dependent bias, and secondly 
by the fact that the SDSS quasar sample only probes the high- 
luminosity tail of the distribution of quasars. 

Several studies have suggested that Civ is the least reli- 
able line to estimate BH mass es (e.g., 'Baskin & Laoi''2005'; 
ISulentic et all 120071: Isiien et a l. 2008; Richards et al. 2011), 
due to a possible non-virial component that contributes to the 
line profile. The model constraints are particularly poor in the 
highest redshift bins. Our rather simplistic model for the re- 
lation between virial masses and true masses may not work 
very well for the problematic Civ estimator A re-assessment 
and possible improvement of the current ver sion of the Civ 
estimator is desirable (e.g., lAssef et ani2011i. Shen et al., in 
preparation). 



We can obtain tighter constraints on the BHMF and model 
parameters if we impose more restrictive prior constraints. 
For instance, if we fix the value of (3, the uncertainties of 
the resulting BHMF and model parameters are substantially 
reduced, while the LF is not changed significantly. This sug- 
gests that better prior knowledge on these parameters can help 
improve these model constraints significantly. 

5. DISCUSSION 

5.1. Evolution of quasar demographics 

It is well know that the numbe r density of bright qu asars 
peaks around redshift '-^ 2-3 (e.g. jRichards et al]l2006l) . The 
most important result in quasar demographics in the past 
decade is the so-called "cosmic downsizing", i.e., the number 
density of less luminous objects peaks at l ower redshift. Ini- 
tially discovered in t he X-ray surveys ( e.g. .ICowie et al.ll2003l : 
ISteffen et alJ I200I lUedaetan I200I iHasinger et all l2005l) . 
this trend is no\y confirmed in the optical band as well (e.g., 
Bongi orno etani2007HCroom et al.ll2009h . 

Here we use the DR7 quasar sample and our model 
LF/BHMF to probe downsizing in quasar demographics. The 
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Fig. 14. — Model BHMF in broad-line quasars. The data points are the binned virial BHMF as in Fig.|3] The magenta shaded region indicates the 68% 
percentile range of the model BHMF for all broad-line quasars, while the green shaded region indicates that for the detectable quasars (e.g., above the flux limit). 
The constraint in a zbin is generally better if there are more quasars in that bin. Note that in general dlogMBH /^log^BH.vir hence the areas underlying the 
green curve and the data points are not necessarily the same. 

flux limit of SDSS quasars is generally not deep enough to 
probe the evolution of the number density of the less luminous 
quasars. However, we have the best constraints on the high- 
luminosity end, and our model LF can extrapolate down to 
fainter luminosities, thus compensating for the shallow depth 
of the SDSS quasar sample. 

Fig.fTSJshows the model LF in all 14 zbins. We only plot 
the portion that is above the flux limit in each individual bin 
for which our model LF reproduces the binned LF well and 
the constraints are the tightest. The width of each stripe in- 
dicates the Icr (68%) range of the model LF. The dashed ver- 
tical line indicates M,(z = 2) = -29.5, beyond which there is 
no binned LF data (see Fig.fTsll. The most noticeable feature 
is that the curvature of the LF changes significantly beyond 
z ^ 3. This is noted in iRichards et al.l (2006) as the flatten- 
ing of the bright-end slope, if a single power-law were fitted 
to the binned LF data. However, it appears to be more of a 
strong evolution of the brea k luminosity than a ch ange in the 
bright-end slope. Recently, Fonta not et al.l (l2007h argue that 
the apparent flattening in the bright-end slope is caused by an 
overestimation of the completeness of SDSS color selection 



at z > 3.5, as was adopted in IRichards et al.l (l2006l) . An im- 
plicit assumption in their argument is that the completeness 
must become lower towards the flux limit, leading to an arti- 
ficial flattening in the bright-end slope. However, we found 
that in our model LF, the break luminosity in the LF gradu- 
ally increases towards higher redshift in a continuous fashion 
(see Fig.fTsTl. and the slope is already significantly flattened at 
z '^ 3.2 for Mi[z = 2] < -27. Thus it is not obvious that this 
flattening in the bright-end slope is caused by the possible dif- 
ferential incompleteness in the SDSS color selection. A more 
careful analysis is needed to resolve this issue. 

Another notable feature is that zbin7 (z = 1.6) shows a 
steeper slope at Mi[z = 2] < -28 than the adjacent zbins. 
This is not caused by the failure of our model, as shown in Fig. 
[T3] A close comparison between our model LF and the binned 
LF is shown in Fig. [16] for zbinG, zbin7 and zbinS. 
The points show the binned LF data, and the lines show the 
model LF from our Bayesian approach; the two are in excel- 
lent agreement. This apparent deviation in zbin7 is likely 
caused by the systematics of emission line /iT-correction. For 
this zbin, the Mgii line contributes the most to /-band. The 
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Fig. 15. — Model LF in the 14 zbins. The width of each stripe indicates 
the 68% range of the model LF. We only show the portion that is beyond 
the flux limit in each zbin and is well constrained by the data. The dashed 
vertical line indi cate s Mj(z = 2) = -29.5, beyond which there is no binned 
LF data (see Fig. [13). The slope for M, [z = 2] < —27 is sig nificantly flattened 
beyond zbinll (1= 3.2). as noted in Rich ards et aljf200Q . It is also notable 
that the slope for zbi n7 is steeper than the adjacent two zbins, which is 
highlighted in Fig. ll6l and discussed in the text. 
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Fig. 16. — Model LF for zbin 6, 7, 8 and comparison to the binned LF 
(open circles with error bars). The solid and dashed lines indicate the median 
and 68% range of the model LF. The binned LF and our model LF are esti- 
mated in completely different ways and are in excellent agreement. The LF 
at z ~ 1.6 has a steeper bright-end slope than at z ~ 1.4 and z ~ 1.8, which 
may reflect systematics in the emission line A'-correction at this redshift (see 
the text for more discussions). 



Mgii line strength relative to the continuum decreases with 
lumin osity (e.g., the B aldwin effect, Baldwin 1977; see fig. 
13 in iShen et al.ll201 1|) , thus using a constant emission line K- 
correction in this bin will lead to more underestimated con- 
tinuum luminosity towards the brighter end. This causes an 
artificial steepening at the bright end of the zbin7 LF(asim- 
ilar effect is seen for zbinl 2 where Civ contributes the most 
to the /-band). This is a second-order effect and therefore we 
do not correct for it in the present work. However, it does 
indicate that as the statistics becomes much better, additional 
systematic effects need to be taken into account in order to get 
unbiased results. We note that zbin7 also shows notable de- 
viations from zbin6 and zbinS in terms of the extrapolated 
LF and constraints on model parameters (see Figs.fTTIandfTSll. 

Finally, as mentioned earlier, the LF (above the SDSS flux 
limit) generally flattens (i.e., the luminosity break) at brighter 
magnitude to wards higher red shift, consistent with earlier 
findings (e.g., ICroom et al.ll2009i) . However, we did not ob- 
serve a significant steepening of t he bright-end slope from 
z = 0.4 to z = 2.65, in tension with ICroom et all (12009). We 
suspect this disc repancy is cau sed by the different methods to 
measure the LF: ICroom et al.l ([2009) fit the LF with a broken 
power-law function, while our LF estimate is non-parametric. 

Fig. \l2\{left) shows the evolution of quasar number density 
as a function of luminosity, using our model LF. Filled cir- 
cles indicate the portion of the LF that is sampled by SDSS 
quasars, and open circles indicates extrapolated model LF. 
Note that at low redshift, the SDSS sample also suffers from 
the bright / = 15 cut, in addition to the faint luminosity cut. 
We only extrapolate the model LF to 3 magnitudes fainter, 
where our model constraints are still reasonably good. With 
the addition of extrapolated data, there is a clear trend that 
the number density of more luminous quasars peak earher. 
Com pared with the results from 2SLAQH-SDSS (ICroom et al.l 
12009). our sample extends to higher redshift and the peaks for 
the bright quasars are well resolved. 

Fig. [it] (ng /if) shows the evolution of quasar number den- 
sity as a function of BH mass. The uncertainties in each 
redshift bin are quite large, reflecting the poorly constrained 
BHMF (e.g.. Fig. \T4\ . However, there is tentative evidence 
that more massive BHs achieve their peak density earlier. This 
is the manifestation of "cosmic downsizing" in terms of BH 
mass. As redshift decreases, active BHs become on average 
less massive and are likely in less massive hosts. The number 
density of Mbh ~ 10^ M0 quasars peaks around z = 2, which is 
consistent with the trend found in ^3.1l based on virial masses 
(although the characteristic BH mass shifted to ^ 3 x lO^M© 
for virial masses due to the luminosity- dependent bias in the 
flux limited sample). iKellv et al] (|201() | ) also found evidence 
for active BH mass downsizing, and IVestergaard & Osmed 
(|2b09) found downsizing at the high-mass end of the active 
BHMF, although the latter did not correct for incompleteness 
and scatter in virial mass estimates. 

We can compare the derive d active BHMF with the local 
dormant BHMF estimated by IShankar et"an (l2009l) . Below 
~ 3 X 10^ Mq (the upper limit that can be measured for the 
local BHMF), the active BHMF at any redshift is always more 
than one order of magnitude lower than the local dormant 
BHMF. This means quasars are by no means a long-lived cos- 
mological population, and we are witnessing different genera- 
tions of quasars at different redshift, which cumulatively (over 
time) build up the local dormant BH population. 

Fig. [TS] further demonstrates the insignificant contribution 
of broad-line quasars to the total BH mass density at all red- 
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Fig. 17. — Downsizing of broad-line quasars. Lefi: the evolution of quasar number density per luminosity interval as a function of luminosity, for our model 
LF. Filled circles are portions of our model LF that ai'e sampled by SDSS quasars, while open circles represent model extrapolations. Error bars stand for the 
68% percentile range from our model LF. The number density of fainter quasars peaks at later time. The glitch at j = 1.6 (zbin7) is discussed in JI5.1I Right: 
the evolution of quasar number density per BH mass interval as a function of BH mass, for our model BHMF (all broad-line quasars). Error bars stand for the 
68% percentile range from our model BHMF. The estimates are quite noisy due to the large uncertainties of our model BHMF. But there is some evidence that 
the number density of lower mass BHs peaks at later time. 
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Fig. 18. — The growth of total BH mass density with cosmic time. The 
black solid line is the accreted BH mass density integrated over the bolomet- 
ric LF determined by Hopkins et al.l |2007|) , assuming a radiative efficiency 
e = 0.1. The gray shaded region has a vertical width of 0.3 dex centered on 
the solid line, which we use as a conservative estimate of the uncertainty of 
the accreted BH mass density due to uncertainties in the bolometric LF and 
radiative efficiency. The black circles are the measured global stellar mass 
density I Perez-Gonzalez et al. 2008), scaled by a factor of 10"^. The ma- 
genta circles are the total BH mass density in broad-line quasars, estimated 
by integrating over our model BHMF. See discussions in the text. 

shifts. The black line shows the accre ted BH mass density 
using the bolometric LF estimated by Hopkins et al.l (|2007|) 
and a radiative efficiency e = 0.1, where the gray shaded re- 
gion is centered on the black line and has a vertical width 
of 0.3 dex, which is a conservative estimate for the uncer- 
tainty in the accreted BH mass density due to uncertainties 



in the bolometric LF and radiative efficiency. The accreted 
BH mass density at z = is consistent with the estimated BH 
mass density using local spheroid-B H scaling relations (e.g., 
iShankar et al.ll2009l:lYu & L"ull2008h . T he black circles are the 
measu red global stellar mass density in lPerez-Gonzalez et al.l 
(|2008'), scaled by a factor of lO""'. Although there are sys- 
tematic uncertainties in both the accreted BH mass density 
and the measured global stellar mass density, the agreement 
between the two is remarkable, and argues strongly for the 
co-evolution between galaxies and BHs. However, the total 
BH mass density in broad-line quasars (magenta circles), es- 
timated by integrating over our model BHMF, is at least one 
order of magnitude less than the total BH mass density at all 
times. 

5.2. The Eddington ratio distribution 

Our model specifies the Eddington ratio distribution (or lu- 
minosity distribution) at fixed BH mass'". Even though the 
constrained parameters have large uncertainties (see Fig.fTTTl. 
it appears that a linear relation between the mean luminos- 
ity and mass is consistent with the data. This means that 
the mean Eddington ratio does not depend on mass strongly. 
However, there is some evidence that the mean Eddington ra- 

'" The Eddington ratio distribution at fixed BH mass is different from that 
at fixed luminosity, in the presence of scatter of the Mbh ~ ^bol relation (cf . 
Fig. [8). This is true no matter which mass is used (virial versus true). For 
the virial mass-based Eddington ratio at fixed luminosity, the average value 
scales as L" '' oc Mgn.vir from the virial relation (assuming a slope of 0.5 in 
the mean R-L relation). However, for the virial mass-based Eddington ratio 
at fixed virial mass, the average value has a much weaker dependence on 
MgHvir. and one also must account for the selection effect of the flux limit. 
These points are well demonstrated in Fig. [8] 
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Fig. 19. — Evolution of the sample-averaged Eddington ratio. The magenta 
points are for all broad-line quasars with true masses; the green points are for 
the detectable quasars with true masses; both are estimated from our model. 
The black points are the observed average Eddington ratio in our quasar sam- 
ple, estimated using virial BH masses. Note that the error bars here indicate 
the uncertainty in the mean value, not the scatter in Eddington ratios. 



tio increases with redshift, which may reflect the evolution 
of the average accretion efficiency or the general decline of 
quasar light curve. The scatter in the Eddington ratio distribu- 
tion is poorly constrained, and no coherent redshift evolution 
is seen. Nevertheless, we found a median value of ^ 0.4 dex 
in the dispersion of Eddington ratios at fixed BH mas s, con- 
sisten t with earlier studies (e.g., Shen et al. 2008; Kell v et all 
l2OT0h . 

The Eddington ratio distribution in any flux-limited sam- 
ple or at fixed luminosity is essentially different from that of 
all broad-line quasars or at fixed BH mass, as emphasized in 
Uhen et al. (2008) and Kelly et al. (2010). It suffers from the 
Malmquist-type bias such that more high-Eddington ratio and 
low mass objects are scattered into the flux-limited sample 
than low-Eddington ratio and high mass objects scattered out 
(i.e., due to a non-zero ct/ and a bottom-heavy BHMF). Thus 
the average Eddington ratio is biased high in our sample. In 
addition, using virial masses can further modify the sample- 
averaged Eddington ratio due to the luminosity-dependent 
bias (i.e., a non-zero /3). 

Fig. [19] summarizes these behaviors. Plotted here are the 
sample-averaged Eddington ratio and its uncertainty" in 14 
zbins under different circumstances. The filled magenta 
squares are for all the broad-line quasars, and the open green 
squares are for those that are detectable above the flux limit; 
both are estimated using our model outputs. The difference 
between the mean Eddington ratio for all quasars and de- 
tectable quasars reflects the scatter in luminosity at fixed true 
mass (ct/) and the shape of the underlying BHMF, and in gen- 
eral the mean Eddington ratio in the detectable sample is bi- 
ased high. The open triangles are for all quasars in our sam- 
ple, where the Eddington ratios are estimated with virial BH 
masses. The difference between the mean Eddington ratio for 
the detectable quasars estimated with true masses and virial 
masses mostly reflects the luminosity-dependent bias (i.e., a 
non-zero /?). Since we have a positive [3 for most zbins, 
the mean Eddington ratio based on virial masses for the de- 
tectable quasars is generally biased low than the true value'^. 

' ' Note that the uncertainty here refers to the uncertainty in the mean Ed- 
dington ratio, not the scatter in Eddington ratio. 

'^ When /3 is small, the difference between the mean Eddington ratio with 



5.3. Model caveats and future perspectives 

Our forward Bayesian modeling stands for a major 
improvement on retrieving information from the mass- 
luminosity plane compared with previous studies based on 
virial BH masses. However, there are still several caveats in 
our approach, which can be improved in future investigations: 

• First and foremost, we assumed that the specific forms 
of single-epoch virial mass estimators adopted in this 
work give unbiased BH mass estimates when averaged 
over luminosity at fixed true mass, i.e., E{mc\m) = m, 
where £(...) stands for the expectation value. While this 
is true for the local RM AGN sample (by calibration), 
it may not hold when extrapolating to the high-redshift 
and high-luminosity population. 

Also, we recognize that different versions of virial es- 
timators (even for the same line) do not necessarily 
yield the same BH masses. If the virial mass esti- 
mators that we used were already biased for objects 
in our sample (i.e., due to an imperfect virial calibra- 
tion, or other generic systematics in the virial tech- 
nique) then the BHMF estimation would also be bi- 
ased. In particular, it is possible that FWHM is a 
worse indicator for the virial velocity than other line 
width mea sures (s uch as line dispersi on, see arguments 
by Peters on et alj|20 04; ColH n et all [2006; Wan g et alj 
l2009u iRafiee & Halli201 la.b.) . and the luminosity bias 
we found here may reflect the systemics with incorrect 
forms of virial estimators. 

Nevertheless, there is currently no consensus on which 
versions of virial estimator are the best. As a proof of 
concept, our model can easily be adapted to use the up- 
dated virial calibrations when they become available. A 
similar study based on line dispersion-based virial BH 
masses is currently under way. 

• Our parameterized model is still in a relatively restric- 
tive form, especially for the the Eddington ratio model 
at fixed true mass. We have tested with rather relaxed 
parameterizations for the joint distribution of luminos- 
ity and true mass, but found that the current data is not 
sufficient to give reasonable constraints if we treat the 
luminosity-dependent bias as a free parameter. This can 
be seen in Fig. [8] that the current data only sample the 
tip of the distribution, and thus it is difficult to con- 
strain parameters for overly flexible models. Deeper 
data is needed to test more flexible models. A fainter 
broa d-line quasar samp le, such as the 2SLAQ sam- 
ple dCroorn et alJ l2009l) or the BOSS quasar sample 
dRoss et al.l201 iV will also provide more stringent con- 
straints on our model parameters and better representa- 
tion of the mass-luminosity plane. 

Nevertheless, the capability of reproducing the ob- 
served distribution is reassuring that a single lognormal 
Eddington ratio model with a mass-dependent mean is 
still a reasonable choice. An alternative model would 
be a power-law Eddington ratio distribution at fixed 

true masses and virial masses for the flux-limited sample may not be discern- 
able, due to the slightly different calculations of the mean Eddington ratio 
in the two cases: the mean true Eddington ratio is calculated using random 
draws from the posterior distribution; the mean virial Eddington ratio is the 
median in the observed distribution. 
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true mass, as suggested in some models fo r the de- 
caying part of the quasar light curv e (e.g., lYu & Lul 
12004 iHopkins et al.ll2005l;IShenl2009l) . However, since 
quasars are continuously forming (especially for fainter 
quasars, whose number density peaks later in time), at 
each epoch we should also witness quasars in their ris- 
ing part of the light curve, unless this period is com- 
pletely obscured. In addition, the broad line region 
may only exist when the BH is accreting at high ac- 
cretion r ate during AGN evolution (e.g., IShen et all 
r2007a; Hopkins etal.ll2009l). i.e., w ith Eddington ratio 
A ~ 10"^- 1 (e.g.. lTrump et alJ201 ll) . and there must be 
fluctuations in the instantaneous Eddington ratio. Thus 
a lognormal model may indeed capture the basic char- 
acteristic of the Eddington ratio distribution at fixe BH 
mass. 

In a follow-up paper (Kelly & Shen, in preparation), we 
will use a more generic prescription to model the joint 
distribution of luminosity and BH mass directly instead 
of using a specific Eddington ratio model as adopted 
here. This will be achieved at the cost of simplifying 
other prescriptions in the model, for instance, fixing the 
value of [3. However, by relaxing the single lognormal 
Eddington ratio model, we can test if a single lognormal 
model is a good approximation. It also allows a more 
proper investigation of the mass-luminosity plane. 

Finally, we note that we assumed a Gaussian distribu- 
tion for the error in the virial mass estimates (which 
may not be adequate), and our model is not robust to 
outliers. These issues may be partly responsible for the 
large uncertainties we get in several zbins. We plan 
to investigate these issues in future work. 

i Since we are only concerned with broad-line objects, 
the active BHMF derived in this paper is a lower- 
limit for the active SMBH population. It is impor- 
tant to include the growth of SMBHs in obscured ac- 
cretion or with accretion mode that does not produce 
broad lines, as traced by populations of active BHs 



selected by various other techniques (e.g.. IStern et al. 



2009; 



2011 



selected oy various otner tecnniques (e.g.. latern et ai 
20051 : iReves etani2008t iLuo et al.ll2008l:lTreister et al 



Hickox et all 120091; lElvis et al.l I2009I;" lYan et al 
; IXue et all 1201 U ICivano et all 1201 ll) . This non- 
broad-line population could contribute as much as one 
half of the total accreted BH mass density, although sig- 
nificant uncertainty remains in determining their rela- 
tive abundance to the broad-line population. 

6. CONCLUSIONS 

This paper represents our first step towards a systematic in- 
vestigation on the demographics of broad-line quasars in the 
mass-luminosity plane and its redshift evolution. We used a 
forward modeling Bayesian framework to model the joint dis- 
tribution in the mass-luminosity plane. With simple model 
prescriptions for the underlying active BHMF and Eddington 
ratio distribution, we were able to fit the observed distribu- 
tion above the sample flux limit and extrapolate below using 
constrained model parameters. We paid particular attention to 
the distinction between virial mass estimates and true masses, 
and corrected for the selection effects of flux limited samples. 
The main conclusions of the paper are the following: 



1 . Virial BH masses are not true masses 03.2.11 ). and their 
explicit dependence on luminosity has implications for 



their conditional probabihty distributions: p{me\m) ^ 
p{me\m,l) ^ p{me\l), and cTvir > cr/,„. We found evi- 
dence that for Mgii and possibly Civ, there is a pos- 
itive luminosity-dependent bias in virial masses, such 
that at fixed true mass, objects with luminosities above 
average have over-estimated virial masses. This is ex- 
pected as there must be uncorrelated scatter between 
line width and luminosity at fixed true mass, due to the 
imperfectness of them being the surrogates for the virial 
velocity and the BLR radius in the virial mas s tech- 
nique . While this bias wa s noted earlier (Shenet al] 
I2008t IShen & KeUyl 1201 Ol) . fliis is for the first time 
that we quantified the level of this bias with the SDSS 
quasar sample. As a consequence, the dispersion in 
virial mass in narrow luminosity bins or flux-limited 
sample s is generally smaller than th e virial mass uncer- 
tainty (IShen et al.l 120081: IShen & K eflv 2010), and the 
average virial BH masses in z > 0.7 SDSS quasars are 
likely overestimated by a factor of a few. 

2. Our model LF is tightly constrained in the regime sam- 
pled by SDSS quasars, and makes reasonable predic- 
tions when extrapolated ^ 3 magnitudes fainter, as 
compared with the LF measured from faint quasar sur- 
veys. Downsizing in LF evolution is clearly seen with 
our model LF. The break luminosity (as in a double 
power-law model) increases towards higher redshift. 
As a consequence, the LF slope for M,[z = 2] < -27 
from a single power-law fit appears shallower at z > 3 
than at lower redshift. 

3. Except for a few zbins, the model active BHMF usu- 
ally cannot be well constrained to within a factor of a 
few. This is mainly caused by the additional free pa- 
rameter on the luminosity-dependent bias, which leads 
to degeneracies with other parameters. The SDSS sam- 
ple only probes the tip of the active SMBH population 
at high redshift, which makes it more difficult to con- 
strain model parameters for the whole population. Nev- 
ertheless, the abundance of the most massive quasars 
is always overestimated using virial masses regardless 
of the luminosity-dependent bias. Within our model 
a turn-over at the low-mass end of the true BHMF is 
needed in order to fit the observed distribution in the 
mass-luminosity plane, but we caution that this result is 
based on our model assumptions and may suffer from 
the poor sampling of low-mass BHs in the SDSS sam- 
ple. Further investigations are needed to test the ro- 
bustness of this feature. This turn-over shifts to larger 
BH masses for the detected population in SDSS due 
to the flux limit. Although with large error bars, we 
found tentative evidence that downsizing also manifests 
itself in the active BHMF, and the number density of 
Mbh ~ 10'* M0 quasars peaks around z ^2. 

The total BH mass density in broad-line quasars is al- 
ways insignificant compared with that in all SMBHs at 
any redshift (Fig.lTsll. 

4. Within our model uncertainties, the Eddington ratio 
distribution at fixed true BH mass has a mean value 
that weakly depends on mass, and an average scat- 
ter of ^ 0.4 dex. The sample-averaged Eddington ra- 
tio for all broad-line quasars ranges between ~ 0.01 
and '-^ 0.3, and increases with redshift (Fig. [T9] l. The 
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sample-averaged Eddington ratio for quasars above the 
SDSS flux limit is Malmquist-biased due to the scatter 
in Eddington ratios at fixed BH mass and the shape of 
the underlying active BHMF. Furthermore, using virial 
masses for the detected quasars tend to reduce the Ed- 
dington ratio again due to the luminosity-dependent 
bias (Fig. [HI). 

5. Our model reproduces the observed distribution in the 
mass-luminosity plane (with virial masses). The flux 
limit, the scatter between true masses and virial masses, 
as well as the luminosity-dependent bias, change the 
distribution in the mass-luminosity plane significantly. 
Thus any features in the mass-luminosity plane based 
on virial masses must be interpreted with caution. 

Our results highlight the need for understanding the system- 
atics in the virial mass technique. All single-epoch virial mass 
estimators currently utilized in the literature are bootstrapped 
from the local RM AGN sample of only several dozen objects. 
These mass estimators have been extensively used for high 
redshift and high luminosity quasars, a regime that is poorly 
sampled by RM AGNs. We have shown that even without 
other systematic issues with the virial technique, the scatter 
and luminosity-dependent bias between virial masses and true 
masses akeady make it difficult to constrain the active BHMF 
accurately. Any additional systematic issues will simply make 
the situation worse. 

Since BH mass is a crucial physical quantity and is re- 
lated to many fundamental processes of SMBHs, refining 
the techniques that weigh the BH is of tremendous impor- 
tance. As the most promising method to measure active 
BH masses, the RM technique (and its extensions of single- 
epoch virial estimators) should be tested and calibrated care- 
fully. Progress is being made in reverberation mapping stud- 



ies of broad -line AGNs and calibrations of virial mass esti- 
mator s (e.g..lKaspi et al.ll2007l: iBentz et aljr2009bt IWoo et alJ 



l2010t lGrahametal.ll2011h . as well as in the statistical de 
scription o f AGN variabil i ty as applied to reverberation map- 
ping (e.g., IZu et alJl20Tlt iBrewer et aDl2011l: iPancoastetaP 
2011). However, to fully understand the BLR kinemat- 
ics/geometry and to construct reliable mass estimators, a sub- 
stantially larger RM sample is needed to sample an unbiased 
parameter space of broad-line objects, as well as a much bet- 
ter understanding of the systematics in the RM technique and 
single-epoch mass estimators. 
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